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ABSTRACT 

This paper describes a series of 3D simulations of shallow inefficient convection 
in the outer layers of the Sun. The computational domain is a closed box containing 
the convection-radiation transition layer, located at the top of the solar convection 
zone. The most salient features of the simulations are that: i)The position of the lower 
boundary can have a major effect on the characteristics of solar surface convection 
(thermal structure, kinetic energy and turbulent pressure). ii)The width of the box 
has only a minor effect on the thermal structure, but a more significant effect on the 
dynamics (rms velocities). iii)Between the surface and a depth of 1 Mm, even though 
the density and pressure increase by an order of magnitude, the vertical correlation 
length of vertical velocity is always close to 600 km. iv) In this region the vertical 
velocity cannot be scaled by the pressure or the density scale height. This casts doubt 
on the applicability of the mixing length theory, not only in the superadiabatic layer, 
but also in the adjacent underlying layers, v) The final statistically steady state is not 
strictly dependent on the initial atmospheric stratification. 
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1 INTRODUCTION 

It is now just possible to perform physically realistic 3D 
simulations of the surface layers of the Sun which take into 
account the complex interaction between radiative and con- 
vective energy transports (Stein & Nordlund 1998, Kim & 
Chan 1998 hereafter denoted SN and KC, respectively, Stein 
& Nordlund 2000). There have also been a number of 2D 
simulations of the surface layers (Steffen et al. 1990 and 
Gadun et al. 2000). To model solar convection realistically, 
requires a realistic equation of state, realistic opacities and 
proper treatment of radiative transfer in the shallow layers. 
SN and KC are the two most freqently cited 3D models of 
this type of stratified convection. As their approaches dif- 
fer considerably, both in numerical methods and in input 
model physics, it is important to find out, what particular 
aspect of the respective simulations, caused their results to 
be different. 

The aim of this paper is to describe solar surface convec- 
tion which not only has the KC realistic physics, but also has 
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a realistic (as is presently possible) geometry. To achieve this 
we had to increase both the depth and width of the original 
KC model, until the side walls and the horizontal boundaries 
had only a minimal effect on the flow. The simulations them- 
selves model a region less than a few thousand kilometers in 
depth, as measured inwards from the visible solar surface. 
In the deep regions of the solar convection zone, the turbu- 
lent velocity is subsonic and the superadiabaticity is close to 
zero. However, within a few hundred kilometers of the solar 
surface, the convective flux starts to decrease. As the total 
flux is constant, the radiative flux must increase to offset the 
drop in the convective flux. This is achieved by a rise in the 
local temperature gradient V. This region of inefficient con- 
vection is called the superadiabatic layer (SAL). In the SAL, 
the superadiabaticity V — V a d is positive and of order unity 
(Demarque, Guenther & Kim 1997, 1999). As the buoy- 
ancy force is large in the SAL, the region is characterised 
by highly turbulent velocities and large relative thermody- 
namic fluctuations. The turbulent velocity also gives rise to 
a significant turbulent pressure (approximately 15 % of the 
gas pressure). This moves out the convection surface, mod- 
ifying the SAL and the stratification. In one-dimensional 
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(ID) models of the solar convection zone based on the mix- 
ing length theory (MLT)(B6hm-Vitense 1958), the veloc- 
ity is set to zero above the convection boundary. However, 
three-dimensional (3D) numerical simulations described in 
Cattaneo et al. (1990), have shown that just above the top 
of the convection layer the turbulent velocities are still high. 

There are several motivations for such simulations 
among stellar physicists. One is to understand the effects 
of turbulence on the structure of the outer solar layers as 
revealed by the observed frequencies of solar p-modes. An- 
other is to explain the excitation mechanism of the p-modes. 
Still another is to derive more realistic surface boundaries 
for stars with convection zones from first physical princi- 
ples, free of the arbitrary assumptions of the mixing length 
theory. And finally, such simulations may be of help in in- 
vestigations of the solar dynamo. 

Over the last few years, the science of helioseismology 
has provided some very precise measurements of the p-mode 
oscillation frequencies (to within one part in a thou sand) 
in the surface layers of the Sun ( Harvey et al. 1996 ). The 
discrepancy between the observed p-mode frequencies and 
those calculated from solar models, is known to be primar- 
ily due to the approximations made in modelling the sur- 
face layers, where turbulent and radiative losses are signifi- 
cant (Balmforth 1992; Guentherl994). In the 3D simulation 
described in KC, the turbulent pressure pushed the con- 
vective boundary radially outwards from its original posi- 
tion (that was computed using the MLT). This situation 
was mimicked in the I D solar model by tweaking the opac - 
ity in the outer layers (Demarque, Guenther & Kim 1999). 
This resulted in improved p-mode frequencies for low and in- 
termediate degrees. However, Demarque et al. also showed 
that the mixing length prescription of Canuto & Mazzitelli 
(1991), which had a completely different SAL structure than 
KC, could produce a similar improvement of the p-mode 
frequencies for the same Z-values. Full details of the differ- 
ent approaches, the contrasting SAL structures and the re- 
sulting p-modes, are described in Demarque et al. (1997, 
1999). Later, Rosenthal et al (1999) used another approach 
to compute the p-mode frequencies. These authors patched 
the mean stratification (horizontal average) of a 3D hydro- 
dynamical simulation, onto a ID MLT envelope model. In 
order to get a smooth ID model, they adjusted the mixing 
length and the amplitude of the turbulent pressure to match 
the 3D simulation. The computed frequencies were closer to 
the observed frequencies than if a standard solar model is 
used, thus showing the importance of including turbulence in 
modeling the outer layers of the Sun. More recently, Li et al. 
(2002) found similar results to Rosenthal et al. by inserting 
the averaged turbulent pressure and turbulent kinetic energy 
directly into the ID stellar model. Their method is applied 
to two of the simulations described in this paper (see section 
4.3). As the turbulent kinetic energy and turbulent pressure 
were very small at the base of the 3D model, they were set 
to zero in regions of the ID stellar model that lay below the 
3D model domain. This required the usual adjustment of 
the mixing length parameter and the helium abundance to 
calibrate the perturbed stellar model. No other adjustable 
parameters were employed. The improvement in the eigen- 
frequencies was found to be primarily due to the inclusion 
of turbulent kinetic energy flux in the ID stellar model. 

After extensive testing, we found that our simulations 



are in good agreement with other numerical studies of the 
surface layers (e.g. Rosenthal et al. 1999, Asplund et al. 
2000). In addition, by incorporating the computed 3D turbu- 
lence into a ID stellar model, we were able to produce solar 
surface eigenfrequencies (p-modes) that were very close to 
the observed frequencies. As our eventual goal is to simulate 
the SAL in stars other than the Sun, it is vital to be sure 
we are modelling the Sun as correctly as possible. 



2 MODELLING REALISTIC SOLAR SURFACE 
CONVECTION 

To model surface layer convection in the Sun as realistically 
as possible, we take the following approach: 

(i) Using a stellar evolution code (YREC e.g. see Guen- 
ther et al. (1992)), we compute a standard stellar model from 
which the initial density p and internal energy e, required 
by the 3D simulations, are derived. From an arbitrary ini- 
tial velocity field v , we then compute p, E(— l/2pv 2 + e), 
pv x , pv y and pv z . These are the dependent variables of the 
governing equations. The horizontal directions are x and y, 
and z is radially outwards. 

(ii) Using the same tables for the equation of state and 
the opacities as in the stellar model, we then compute 
the pressure P(p,e), temperature T(p, e), Rosseland mean 
opacity K,(p,e), specific heat capacity at constant pressure 
c p (p, e), adiabatic gradient V a d(p, e) and some thermody- 
namic derivatives. 

(iii) The radiation flux is then computed using the diffu- 
sion approximation in the optically thick regions and the 3D 
Eddington approximation in the optically thin layers. 

(iv) We then integrate the Navier-Stokes equations over 
one time step to compute a new set of dependent variables 
and return to (ii). 



2.1 Realistic initial conditions: steps (i) and (ii) 

The solar model uses the same realistic physics as described 
in Guenther & Demarque (1997). In particular, the low tem- 
perature opacities of Alexander & Ferguson (1994) and the 
OPAL opacities and equation of state were used (tglesias 
*k Rogers 1996). Hydrogen and helium ionisation, and the 
diffusion of both helium and heavy elements are included. 

In the initial model, the atmospheric layers, which are 
in radiative equilibrium, are assumed to be gray. Deeper 
in, in the convectively unstable region, the thermal struc- 
ture is described by the MLT, which prescribes the temper- 
ature gradient V. In the original KC simulation, the Ed- 
dington approximation T(r) relation was used in the at- 
mosphere. In this case, the values of the parameters, X, Z 
and a, in the calibrated Standard Solar Model (SSM) are 
(X,Z,a) = (0.7385,0.0181,2.02), where X and Y are the 
hydrogen and helium abundances by mass, and a is the ra- 
tio of mixing length to pressure scale height in the convection 
zone, required to match precisely the solar radius. It should 
be mentioned that although the values of X and Y in the 
calibrated SSM depend little on the treatment of the sur- 
face layers, the value of a is sensitive to the choice of T(t) 
relation in the atmosphere (Guenther et al. 1992). 

We note that for simulations discussed in this paper 



Solar surface simulations 3 



(those listed in the append ix), the empirical Krish na- Swamy 
T(t) relation for the Sun (Krishna Swamy (1966)) was used 



instead of the Eddington approximation. The same model 
as KC, but constructed with the Krishna-Swamy T(r) re- 
lation, yields (X,Z,a) = (0.7424,0.01706,2.1319) (this is 
model KC2 in the appendix). Because of the inclusion of he- 
lium and heavy element diffusion during the evolution, the 
initial X and Z (denoted by X0,Z0) were slightly different. 
For KC they were (X0,Z0) = (0.7066, 0.0201), and for KC2 
they were (X0,Z0)=(0.710,0.019). The difference in initial 
atmospheric structures between simulations KC and KC2 
provides us with the opportunity to verify that that the fi- 
nal state in the simulations is not strictly dependent on the 
choice of T(r) relation in the initial model. This important 
test is discussed in the appendix. 

For each timestep of the numerical integration we need 
to solve the complex equation of state (step (ii)). This adds 
considerable computational time compared to an ideal sim- 
ulation, such as with a perfect gas. The 3D hydrodynamical 
simulations use identical opacities and equation of state as 
used in the ID reference SSM, which served as the initial 
model. The side boundaries are periodic, while the top and 
bottom boundaries are stress free. A constant heat flux flows 
through the base, and the top is a perfect conductor. To en- 
sure that mass, momentum and energy are fully conserved, 
we use impenetrable (closed) top and bottom boundaries. 

A particular model is specified by g (the surface gravity) 
and T e B (the effective temperature). Aside from the viscosity 
coefficients there are no other free parameters. 



2.2 Radiative transfer: step (iii) 

In the SAL, the photon mean free path may not be small 
enough to use the diffusion approximation. Consequently 
one is forced to either solve the full radiative transport 
equation or consider the three-dimensional Eddington ap- 
proximation (Unno & Spiegel 1966), which is a higher or- 
der approximation than the diffusion approximation, and 
is valid in the optically thin regions. Computationally, to 
solve the full radiative transport equations is formidable, 
and only a few ray directions are currently used in this ap- 
proach (SN). In our simulations, we have chosen to use the 
three-dimensional Eddington approximation. 

In the deeper part of the domain (r > 10 4 ), we use the 
diffusion approximation, 
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where re is the Rosseland mean opacity, a is the Boltzmann 
constant and c is the speed of light. In the shallow region 
Q ra d is computed as 
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where the mean intensity J was computed by using the gen- 
eralized three-dimensional Eddington approximation (Unno 
& Spiegel), 



V ■ | — - VJ ) - npj + k P B = 0, 



(3) 



where B is the Planck function. This formulation is exact 
for isotropic radiation in a gray atmosphere, and without 



requiring local thermodynamic equilibrium, the Eddington 
approximation describes the optically thick and thin regions 
exactly (Rutten 1995). To study spectral line profiles and 
spectral energy distribution, requires frequency dependent 
radiative transfer. However, in a study of the SAL, to serve 
as a surface boundary condition for stellar models or for 
comparison with results of helioseismology, we found a gray 
atmosphere to be adequate. 



2.3 Hydrodynamics: step (iv) 

For deep (v « c s , where v is the flow velocity and Cs is the 
isothermal sound speed) and efficient (V — V a d just above 
zero) convection, Chan & Sofia (1989) showed that the MLT 
is a very good approximation to the real situation. However, 
in the SAL, both v/ca and V — V a d, can be of order unity. In 
this case, the MLT is unlikely to apply. The validity of this 
argument is confirmed by helioseismology. The run of the 
sound speed in the SAL derived by inversion of the helioseis- 
mic data does in fact disagree with the MLT model (Basu & 
Antia 1997). In such an environment, the governing hydro- 
dynamic equations are the fully compressible Navier-Stokes 
equations (see for example Kim et al. 1995). 



dp/dt 
dpw/dt 
dE/dt 



-V ■ pv 



—V ■ pvv 

-V • [(E + P)v - v 
f pv ■ g + Q rad 
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where E = e + pv /2 is the total energy density and p, v, P, e 
and g, are the density, velocity, pressure, specific internal 
energy and acceleration due to gravity, respectively. Q ra d is 
the energy transferred by radiation (see previous section) 
and f is the diffusive flux. Ignoring the coefficient of bulk 
viscosity, the viscous stress tensor for a Newtonian fluid is 
Ey = p(dvi/dxj + dvj/dxi) — 2p/3(V • v)5ij, where p is the 
dynamic viscosity and <5y is the delta function. 

In fully developed turbulence the ratio of the length of 
the largest eddy to the dissipation length is Re 3 / 4 , where Re 
is the Reynolds number, (Landau & Lifshitz 1987). In the 
Sun, Re is the order of 10 12 which means about 10 9 scales 
per dimension. A 3D direct numerical simulation of the Sun 
would thus require about 10 27 grid points ! 

The Large Eddy Simulation (LES) approach assumes 
that the small scales are independent of the resolved scales 
(large eddies) and can be parameterised as a diffusion pro- 
cess. In this case p is an eddy viscosity defined in terms of 
the resolved velocity (Smagorinsky 1963), 



p = p(c M A) (2a : ct) 



1/2 



(7) 



The colon inside the brackets denotes tensor contraction of 
the rate of strain tensor Oij = (Vi«j + VjVi)/2. The sub- 
grid scale (SGS) eddy coefficient c^, is set to 0.2, the value 
for incompressible turbulence, and A = (A x Aj,) 1 '' 2 A z is an 
estimate of the local mesh size. To handle shocks, p is mul- 
tiplied by 1 + C ■ (V • v) 2 , where the constant C is made as 
small as possible, while still maintaining numerical stability. 
As p is dependent on the velocity divergence, any large ve- 
locity gradients are smoothed out by the increased viscosity. 
If a shock occurs it is not resolved, but smeared out by a 
local increase in viscosity. 
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The diffusive flux f = — (/x/Pr)TVS r , where the horizon- 
tal mean of the entropy gradient X7S < i.e. the convection 
zone, and f = — (/xc p /Pr)VT where the horizontal mean of 
V5 > i.e. the radiation zone (the Prandtl number Pr is 
defined below) . In the convection layer the SGS diffusive flux 
tends to smooth out entropy fluctuations and make the layer 
close to adiabatic (in analogy with turbulent mixing). The 
change in the form of the diffusive flux above the convec- 
tion boundary is necessary because the SGS's should con- 
tinue to transport heat radially outwards though the top of 
the box. Away from the horizontal boundaries, f is close to 
zero. At the base, f is equal to crT^j, where a is the Stefan- 
Boltzmann constant. The Prandtl number Pr = v/k, where 
v is the kinematic viscosity and k is the thermal diffusivity. 
In the simulations Pr = 1/3. Due the inclusion of radiative 
energy transport the effective Pr is actually much smaller 
and not constant. 



3 NUMERICAL INTEGRATION: OBTAINING 
ACCURATE STATISTICS 

To simulate the highly stratified SAL of the Sun we need to 
relax the initial layer and then compute accurate statistics. 
The former requires a long computation, while the latter 
a small timestep. In compressible hydrodynamics however, 
with an explicit numerical method, the timestep must be 
less than the time for a sound wave to traverse two adjacent 
grid points. This is known as the Courant-Friedrichs-Levy 
(CFL) stability criterion. Because of these considerations 
the simulations were done in two stages. 

Firstly, using an implicit code in which the time step 
is restricted by the flow speed rather than the sound speed, 
the initial hydrostatic layer was allowed to adjust its ther- 
modynamical structure until it was close to hydrodynamic 
equilibrium. This thermal adjustment phase took at least 5 
hours of solar surface convection time. 

Secondly, using a second order accurate explicit code, 
quantities were averaged over a time that was long enough 
for the averages to be independent of the integration time. 
As the explicit timestep is about five times smaller than the 
implicit timestep, prior to statistical averaging, the code was 
run for a few thousand timesteps. This allowed the simula- 
tion to adjust to the new timestep. The statistical conver- 
gence took at least an hour of solar surface convection time. 



3.1 Thermal adjustment 

The implicit code was the Alternating Direction Implicit 
Method on a Staggered grid (ADISM) developed by Chan 
& Wolff (1982). This code was used to relax the fluid to a 
self consistent thermal equilibrium. 

The entire layer was assumed to be relaxed when 

• the energy flux leaving the top of the box was within 5 
% of the input flux at the base; 

• the horizontally averaged vertical mass flux was less 
than 10 -4 g/cm 2 /s at every vertical level; 

• the overall thermal structure did not change much over 
time; 

• and the maximum velocity in the box was roughly con- 
stant. 



These criteria must be satisfied, before any useful statistical 
data can be gathered. 



3.2 Statistical convergence 

A second order explicit method (Adams-Bashforth time in- 
tegration) gathered the statistics of the time averaged state. 
This code is much more accurate than ADISM, for example 
the mean energy flux leaving the top of the box was within 
1 % of the input flux at the base. On a 667 MHz Alpha pro- 
cessor, each integration step on a 80 x 80 x 80 grid, required 
about 5 seconds of CPU time. 

The time required for statistical convergence depends 
on the particular quantity being averaged. Conserved quan- 
tities converge very fast. For example, the horizontally av- 
eraged vertical mass flux was less than 10 -5 g/cm 2 /s after 
5 minutes of solar time integration. On the other hand, rms 
velocities converged in a few eddy turnover times (at least 30 
minutes of solar time). While second order turbulent quanti- 
ties, such as the horizontal Reynolds stress, took even longer 
to converge (at least 80 minutes of solar time). 



3.2.1 Some statistical definitions 

In a turbulent fluid a quantity q can be split into a mean 
and a fluctuating part, 



q = q(z) +q'(x,y,z,t). 



(8) 



The overbar represents a combined horizontal and temporal 
average, i.e. 



?(*) 



1 



ti—t\ J V (L X Ly) 



qdxdy I dt 



(9) 



ii, is a time after the system has reached a self-consistent 
thermal equilibrium (the thermal adjustment time). L x and 
L y are the horizontal widths of the box in the x and y direc- 
tion respectively. The time required for statistical conver- 
gence is ti—t\. 

The rms value of a quantity q is defined as 

q" =¥-?, (10) 

while the correlation coefficient of two quantities gi and qi , 
is defined as 



C[q[q 2 ] = 



qi q2 — qi qi 



(ID 



As the simulations have periodic side boundaries, sym- 
metry requires that 

. C[v' x v'y] = 
// // 

• V x =Vy 

The run of v' x ' and after 80 minutes of time inte- 
gration is shown in Fig. nl The closeness of the horizontal 
velocities confirms that the simulation is close to statistical 
convergence. 

By examining many different simulations we found that 
the run of C v' x v' y ] is generally a much stricter test of conver- 
gence. Fig. £ shows C[i4t>y] measured for 4 different integra- 
tion times. Even after 80 minutes, C[i4^] is still about 0.1 
near the bottom. Convergence is faster in the upper layers 
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because when the convecting fluid elements move up through 
the stratification, their rapid expansion smoothes out small 
scale fluctuations. 



4 MAIN RESULTS 

4.1 Effect of domain size 

After substantial numerical testing which is described in de- 
tail in the appendix, we found that : 

(i) An undesirable effect of the impenetrable horizontal 
boundary at the bottom was too speed up the overall flow. 
This produced an artificially high convective flux. Just be- 
low the surface the radiative flux is a significant fraction of 
the total flux. As the total energy flux is fixed, to accom- 
modate the increased convective flux, the radiative flux had 
to reduce. This was achieved by an (unphysical) drop in the 
temperature gradient V in the SAL region. 

To avoid this, the lower boundary has to be positioned 
far enough away from the surface, so that the velocity at 
the base is small and uncorrelated from the motions near 
the surface. 

(ii) If the width of the box was too small then the tur- 
bulent kinetic energy of the granules was artificially small. 
This is because the movement of the larger granules was 
restricted by the walls of the box. 

To avoid this, the aspect ratio was doubled until the tur- 
bulent kinetic energy was unaffected by a further increase in 
aspect ratio. For the Sun this required a width of 2.7 Mm. 



4.2 Comparison with previous 3D numerical 
simulations 

4-2.1 Vertical velocity 

Using the SN code, Asplund et al. (2000) computed a series 
of simulations of solar surface convection in a domain with 
a depth of about 4 Mm and width of 6 Mm. The best res- 
olution in Asplund et al. was 200 2 x 82. The rms vertical 
velocity and the mean velocity are shown in Fig.s 3 and 4 
in their paper. Corresponding velocity plots from our best 
model (model C in the appendix) are shown in Fig. |S| Note 
that we have radially upwards as the positive direction, the 
opposite to the Asplund paper. Despite all the differences 
between the SN and KC approaches, away from the bound- 
aries, the rms and mean vertical velocities in our best model 
are very similar to those in the best Asplund et al. simula- 
tion. The most noticeable differences occur at the top and 
bottom because in our simulations the vertical velocity is 
forced dropped to zero. In Asplund et al. the transmitting 
boundaries allowed the velocity to decrease more gently. 



4- 2. 2 Superadiabaticity 

The superadiabaticty V — V a( j for the MLT, and models KC2 
and C, are plotted in Fig. ^. The abscissa is the radius of the 
simulation divided by the radius of the Sun. For model C, 
the superadiabaticity has a maximum of about 0.6, which is 
close to the value given by the SN code (see Fig. 3 in Rosen- 
thal et al.). The position of the top of the convection layer 
(as determined by the Schwartzchild criterion) was pushed 



out further in KC2 than in C. This is because of the higher 
turbulent pressure. 

When the SAL was moved outwards the convective ef- 
ficiency was reduced and radiation was forced to carry more 
of the total flux. This resulted in an increase in the height of 
the SAL in C (triple dot dashed line) , compared to the MLT 
(crosses). However, in KC2 (dashed line) the convection was 
speeded up by the lower boundary and thus is (incorrectly) 
more efficient than the MLT. This resulted in a drop in the 
height of the peak of the SAL compared to the MLT. 

4.3 Comparison with observational results: 
p-mode oscillation frequencies 

4-3.1 Implementation of 3D turbulence into ID stellar 
models 

By analogy with the work of Lydon & Sofia (1995) on mag- 
netic effects, the 3D turbulence is parameterised in terms of 
two quantities (Li et al. 2002), the turbulent kinetic energy 
per unit mass, 

X=\v" 2 (12) 
where v" 2 — v" 2 + v' y ' 2 + v 1 ! 2 is a dimensional quantity, and 



an anisotropy parameter, 
7 = 1 + 2K'A/') 2 . 



(13) 



The z direction is parallel to the radial direction. 

The turbulence is included by incorporating 2 new vari- 
ables 7 and x into the ID stellar model. To understand how 
this is works, consider a perfect gas in which the ratio of the 
specific heats 7 = c p /c v , the internal energy e = c v T and 
the gas pressure P — pRT. The quantities c v and R are for 
unit mass of a gas and T is the temperature. The previous 
three equations can be expressed as, 



7 = 1 + P/(pe) 



(14) 



If we replaced gas quantities by turbulent quantities, 
i.e. P by P turb (= pv z " 2 ) and e by x and rearrange, then we 
would get: 

Pturb = (7 " 1)PX- (15) 

where 7 is defined in terms of turbulent quantities (i.e. rms 
velocities). Including \ and 7 is equivalent to including x 
and Pturb in the ID stellar model. The two variables 7 and x 
are included in the mathematical system in a self consistent 
manner (see Li et al. for the full mathematical treatment) 

For example, the equation of state for the ID model 
becomes, 

p = p(P T ,T,x,7) 

where P T = -P gas + fkd + -Pturb- 

The continuity equation and the equation of transport of en- 
ergy by radiation remain the same regardless of turbulence. 
In terms of 7 and x> the hydrostatic equilibrium is, 

dP T _ GM r 2(7 - l) x . . 

dM r 47rr 4 47rr 3 ■ 1 ' 

where Mr, G and r have their usual standard stellar model 
(SSM) definitions. The energy conservation equation, 

dL r dS , . 

8M= e - T lt> (17) 
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is also affected by the inclusion of x, as 

TdS = dU + (Pt - Pturb)d(l/p) + d X , (18) 

where d\ represents the work done by turbulence. The quan- 
tities L r and e have their usual SSM definitions. 

The convective flux in the ID model F conv = 
pTDS V conv , now includes the turbulent kinetic energy flux 
xV'conv (without turbulence \ — 0, so d% = x)- V'conv is the 
MLT convective velocity and DS is the entropy excess. Equa- 
tion |l8| is the 'non-standard form' of the 2nd law as described 
in Lydon & Sofia. As turbulence is in general anisotropic it 
is wrong to treat the work done by turbulence as — -Pturb dV. 
This is essentially because Pturb is a tensor and P gas is a 
scalar. 

The equation of energy transport by convection, does 
not change in form, but V is different from that without 
turbulence. The equations that govern envelope integrations 
are changed accordingly (see Li et al.). 

4-3.2 Solar p-mode oscillation frequencies 

To investigate the effect of turbulence on the solar p-mode 
oscillation frequencies, the runs of \ and 7 for models KC2 
or C were incorporated into the ID stellar model using the 
Li et al method. For each run, the p-mode frequencies for 
I — 0, 1, 2, 3, 4, 10. ..100, were computed using Guenther's 
pulsation code, (1994), under the adiabatic approximation. 
When comparing the computed and observed eigenfrequen- 
cies, one should really include the proper modeling of ra- 
diative gains and losses by computing the non-adiabatic fre- 
quencies (Guenther 1994). However, to help isolate the ef- 
fects of turbulence on the p-mode frequencies, we computed 
the adiabatic frequencies. The difference between observed 
and computed adiabatic p-mode frequencies for the standard 
solar model is shown in Fig. |j| The frequency difference is 
scaled by the mode mass Q„i (e.g., Christensen-Dalsgaard 
& Berthomieu 1991). 

The Q n i weighting attempts to correct the p-mode fre- 
quency differences, by removing the dependency on mode 
inertia. The mode differences with high mode inertia get 
more weight than those with lower inertia. This tends to 
give greater weight to the low I, deep penetrating modes 
and less significance to the very high I shallow modes. The 
result is the removal of the I dependence in the frequency 
differences due to perturbations in the near surface regions. 
We are left with only n dependence, which is equivalent to 
frequency dependence. The Q n i weighting enables us to see 
that the discrepancy between the observed and the com- 
puted adiabatic frequencies is worse at the high frequencies, 
thus pointing to a problem in the surface layers. Without 
the Q n i weighting, the mismatch would appear to varying 
degree, at all frequencies and therefore it would not be as 
easy to claim that the outer layers are responsible for the 
mismatch. 

In the statistically steady state (after the flow is ther- 
mally relaxed), the turbulent kinetic energy flux should be 
proportional to the energy input rate at the base. The tur- 
bulent kinetic energy \ should be independent of the geom- 
etry of the computational domain. If the box is too shal- 
low or the width too narrow, then \ can be wrong. To 
illustrate the effect of the geometry on the frequencies, 7 
and x from the KC2 model was used to compute a stellar 



model. The computed adiabatic frequencies with turbulence 
derived from KC2 are shown in Fig. ^j. The resulting fre- 
quencies are much worse than the standard solar model. 
However, if we use model C instead of KC2, then the de- 
rived frequencies were improved considerably (Fig. ^). The 
difference between the computed and the observed frequen- 
cies is less than 5 p Hz. The inclusion of turbulent kinetic 
energy in ID models in addition to turbulent pressure, has 
a similar effect on the adiabatic frequencies to the inclu- 
sion of radiative losses/gains (computation of non-adiabatic 
frequencies). Both reduce the difference between the adia- 
batic frequencies and the observed frequencies by an order 
of magnitude. This mean that the turbulent correction to 
the p- modes is as important as radiative losses/gains. 

4.4 Characteristics of the solar granules 

4-4-1 Size 

Fig. |^ shows a set of contours of the instantaneous vertical 
velocity. From top to bottom the frames depict depths of 
-0.14 Mm, 0.03 Mm, 0.2 Mm and 1.0 Mm. The depth was 
measured positively inward from the visible solar surface. 
The contours themselves have been derived from a simula- 
tion of 80 3 grid points, in a domain with a horizontal area 
of 3.75 Mm x 3.75 Mm, and a depth of 2.5 Mm. At this 
instant in time there appear to be between 3 and 4 gran- 
ules near the surface. The thick black lines represent the 
strong downflows which occur at the sides of the granules. 
The lighter regions denote upflowing fluid or weak down- 
flows. The granular pattern seems to persist over the first 
three frames, suggesting that the granule structure is corre- 
lated over most of the SAL (between Mm and 0.25 Mm). 
By the fourth frame there is not much sign of granulation. 
As the last contour is still more than 2 pressure scale heights 
above the base, the break up of the granules is probably not 
due to the impenetrable bottom boundary. In other words 
the bottom boundary is far enough away from the surface. 
Furthermore, as the shape of the granules does not appear 
to be influenced by the side walls, the width of the box is 
also big enough. 

The thin closely spaced dark vertical/slanted parallel 
lines that appear in the first contour, are signs of two grid 
waves. This is an indication of insufficient damping and is 
only seen at the very top of the box. It is due to the low 
SGS viscosity, which is proportional to density ((equation 
Q). This is is an undesirable numerical effect. But as the os- 
cillations occur only at the very top of the box, where the 
density is very small, their effect on the convection below 
is minor. Before starting the time integration, we were able 
to damp out such oscillations, by slightly increasing the vis- 
cosity near the top. As these contours only represent one 
instant in time, such a picture can only provide a limited 
idea of the nature of the solar granulation. If the same con- 
tours are computed at a later time, a completely different 
picture would be seen. For example the third frame of figure 
22 is the same simulation, but the contour was computed 2 
solar minutes later. 

Clearly any useful data can only be derived from long 
time averaged statistics. One characteristic vertical length 
scale associated with convective turbulence is the half- width 
of the 2-point vertical velocity correlation C^vi]. If we as- 
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sume convective eddies to have an aspect ratio of unity, then 
this length scale gives us an idea of the size of the eddies in 
a turbulent fluid. Using simulation D (described in the ap- 
pendix), we computed Cfw^w^] at a series of depths from the 
top of the box inwards. The results are shown in Fig. |^. 
From about Mm to 1 Mm the eddy size remains close to 
600 km. From the begining to the end of the plateau, both 
density and pressure have increased by an order of magni- 
tude. Over this range the eddy does not seem to be affected 
by stratification. However, over the next Mm the eddy size 
increases to about 1100 km. After that the lower boundary 
starts to influence the eddies. The bends at the far left and 
far right of the plot signal the approach of the upper and 
lower boundaries. 

The velocity contours and the velocity correlation 
length both suggest that the granules have some coherent 
vertical structure. The contours show intergranule lanes that 
do not shift much between and 0.2 Mm, i.e. the gran- 
ules are like cylinders. The half-width suggests a vertical 
size of about 600 km, between 1 Mm and the surface. This 
length scale does not seem to be affected by the stratifica- 
tion. The strong vertical coherence implies that the convec- 
tion is dominated by strong downflows that originate close 
to the surface. SN showed that solar granulation is primarily 
driven by radiative cooling at the surface. When ascending 
fluid approaches the surface, it looses heat so rapidly that 
very strong downflows are created. The downflows are not 
deflected by the surrounding fluid until they are about 1 
Mm below the surface, by which point they have weakened 
enough to be deflected/broken up by the surrounding fluid 
motion. 

4-4-2 Heat transport 

In Chan & Sofia (1987,1989) the 2-point vertical velocity 
and temperature correlations were found to scale with pres- 
sure scale height rather than density scale height. This study 
was for deep efficient convection. In shallower layers Kim 
et al. (1995) found that while the vertical velocity correla- 
tion scaled with both density and pressure scale height, the 
temperature could not be scaled with either. The difference 
from Chan & Sofia's results, was caused by the inclusion of 
the coupling of the partial ionisation with convection (they 
were treated separately by Chan & Sofia). As a fluid parcel 
moves upwards through regions of decreasing ionisation, it 
liberates ionisation energy which increases the buoyancy of 
the fluid parcel. However, because radiation was modelled 
by the diffusion approximation (which is known to break 
down at some point in the SAL) , the Kim et al. simulations 
could only include the lower half of the SAL. 

The present models include all of the SAL. Figs, juj and 
|ll|show C[v' z v' z ] and C[S'S'] at depths of 0.06 Mm, 0.2 Mm, 
0.3 Mm, 0.47 Mm and 1 Mm, denoted by solid, small-dashed, 
dot-dashed, triple-dot dashed and long-dashed lines, respec- 
tively. The curves have been centered about their respec- 
tive maxima. As the plots do not coincide, neither quantity 
scales with the pressure scale height. Furthermore, if den- 
sity is used rather than pressure, there was no noticeable 
improvement (not shown). 

The width of the C[v' z v' z ] distribution decreases between 
0.06 Mm and 1 Mm. This is because the width of the geo- 
metric distribution is roughly constant (i.e. the half-width is 



always about 600 km, as shown in the previous figure), but 
the scale height increases inwards. The entropy correlation 
looks very different. For AlnP < the entropy correlation 
for the three most shallow depths drops rapidly, while for 
A In P > it shows a smoother decay. This implies that as- 
cending parcels near the surface, lose heat very quickly, while 
descending parcels maintain their heat content appreciably 
longer. Over the last Mm before the surface, the eddies may 
have about the same diameter, but they transport less and 
less heat as the surface is approached. 

This behaviour can be understood by examining the 
correlation C[v' z S'] for upflows and downflows. Fig. [l2| shows 
that the granules are more efficient at transporting entropy 
downwards than upwards. At the peak of the SAL, Cfu^S 1 '] 
is 0.9 for downflows and 0.6 for upflows. This explains the 
one-sidedness of C[S'S']. Consider a fluid parcel that is ini- 
tially just below the surface (say at a depth of 0.1 Mm). 
If that parcel moves downwards it will maintain its origi- 
nal entropy for about 1 PSH, whereas if the parcel moves 
upwards it it will loose its original entropy in 0.5 PSH. At 
greater depths C^S 1 '] is similar in the upflows and down- 
flows, hence C[S'S'] is more symmetric deeper down. 

4-4-3 A variable mixing length ? 

Over the last 1000 km before the surface, the eddies seem to 
have a nearly constant diameter of about 600 km. This would 
seem to imply a mixing length ratio (the correlation length 
divided by the pressure scale height) which increases towards 
the surface. However, as the surface is approached, radiation 
plays a greater and greater role in heat transport, so the 
mixing length ratio should approach zero at the convection 
surface. This paradox can be partly resolved by noticing that 
as the surface is approached, the eddies become less and less 
efficient at transporting heat upwards. This is shown in the 
plot of C[i^S'] for the upflows. A more reasonable candidate 
for the mixing length might be the product of C[v' z S'] and 
the velocity correlation half width. 

5 SUMMARY AND IMPLICATIONS FOR 
SOLAR MODELS 

As we intend to use our code to model the SAL in other stars, 
this paper is an important benchmark for future studies. 
From these simulations of the Sun, we found that: 

(i) The impenetrable lower boundary needs to be far 
enough away from the SAL so that by the time the fluid 
reaches the boundary, the velocities have become both weak 
and uncorrelated from velocities in the SAL. If the lower 
boundary is too close to the SAL, the kinetic energy will be 
overestimated. 

(ii) The horizontal cross section needs to be large enough 
so that the side-walls do not restrict the movement of the 
granules. If the box width is too small then the kinetic energy 
will be underestimated. 

(iii) There is a region close to the surface in which the 
vertical correlation length remains constant, even though 
the density and pressure vary by an order of magnitude. 

(iv) In that region the mixing length theory, which as- 
sumes the correlation length to be a constant multiple of 
the pressure scale height, will not work. 
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(v) The final equilibrium state is not strictly dependent 
on the initial model atmosphere (see appendix, section 3.1). 



While these results have been found in a simulation 
of the Sun, it is reasonable to assume that similar criteria 
would apply to the convection-radiation transition layers in 
other stars. The effect of the boundaries should certainly be 
considered when simulating other convection-radiation lay- 
ers. The correlation length (half-width) seems to be a very 
robust feature of the solar granules and could easily be mea- 
sured in other computations. For example, in a recent pro- 
ceedings (Robinson et al 2002) we describe the application 
of this model to the Sun at the sub-giant and the start of the 
red giant branch. Preliminary results indicate that, near the 
surface of each convection zone, the ratio of the half-width 
to the stellar radius depends directly on the surface gravity. 
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ment of Astronomy, Yonsei University grant (2001-1-0134). 
DBG's research is supported in part by a grant from NSERC 
of Canada. We also would like to acknowledge computer sup- 
port from the Centre for High Performance Computing at 
Saint Mary's University in Canada. Finally, we thank H.-G. 
Ludwig for helpful comments. 



6 APPENDIX : NUMERICAL TESTS 

Numerical simulations of solar granulation need to be ro- 
bust. Part of this means that the walls of the computational 
box should not control the radiative hydrodynamics of the 
interior (e.g. the size of the granules, the structure of the 
SAL, the turbulent pressure, etc.). Nor should any small 
changes in the initial conditions affect the final equilibrium. 
The following simulations were designed to minimise these 
uncertainties. 

The features of the individual simulations are sum- 
marised in Table |l|. In this table the geometric width in 
both the x and y directions is given in Mm, while the depth 
is given in Mm and as the number of pressure scale heights 
(PSH). Models A and B were created by vertically extend- 
ing model KC2, while C, D and E were made by periodically 
extending model B. After periodic extension, the horizontal 
resolution in models D and E was halved. Each time a new 
model was created, it was allowed to return to thermal equi- 
librium before the statistical integration was begun. 

For each model we computed three non-dimensional sta- 
tistical quantities: 

(i) The turbulent pressure divided by the mean gas pres- 
sure, 



Pturb = pv" 2 /P 



(19) 



where the '*' denotes a non-dimensional quantity. 

(ii) The turbulent kinetic energy per unit mass divided 
by the isothermal sound speed squared, 

(20) 



X ~2"c7 



where v" 2 = v" 2 + v'' 2 + v" 2 and c s = 



P I p is the isother- 



(iii) The superadiabaticity 

V- Vad 



dlnT 



(21) 



mal sound speed. 



dlnP " u 

where V a d is computed using the OPAL equation of state 

From now on we will drop the '*' on Pturb and \, for conve- 
nience. All quantities in this section are non-dimensional. 

6.1 The influence of varying the initial conditions 

Here we compare the original KC simulation with simula- 
tion KC2 listed in Table 1. Both simulations have the same 
input physics (opacity tables, equation of state) . As noted in 
section 2.1, the only difference between KC2 and KC, is that 
the ID initial models were based on slightly different model 
stellar atmospheres [i.e. the Eddington T(r) relation in KC, 
vs. the empirical Krishna-Swamy T(t) relation in KC2]. As 
the Krishna-Swamy relation is from observations of the Sun, 
it was closer to the final state of the solar simulation. 

Fig. [l] shows the initial hydrostatic (based on the MLT) 
and relaxed hydrodynamic runs of Log P vs Log T for the 
two simulations. The dotted and solid lines are the initial 
plots for KC and KC2, respectively. The horizontal and 
temporal averages of the relaxed states, are denoted by dia- 
monds and crosses for KC and KC2, respectively. To demon- 
strate statistical convergence, the data of model KC2 was 
initially averaged over about 600 seconds of solar surface 
convection, whereas in KC, the averaging time was 2240 
seconds. After a sufficient amount of time, KC2 (crosses) 
converged to KC (diamonds). The MLT provides the ini- 
tial thermal structure, but the hydrodynamical turbulence 
shapes the final structure. Provided the initial conditions 
are not too different, the hydrodynamics converges to the 
same equilibrium state. 

This is an important preliminary result. Even if the ini- 
tial atmospheres are slightly different, the hydrodynamics 
dictates the final equilibrium. The thermodynamics of the 
resultant 3D simulation does not depend on the precise de- 
tails of the MLT. The final equilibrium is not determined 
by the exact value of the mixing length ratio in the stel- 
lar model, on which the initial state of the simulation is 
based. The choice between the two initial conditions in the 
atmosphere, does not affect the turbulent pressure, turbu- 
lent kinetic energy, or the superadiabatic temperature gradi- 
ent. These quantities will depend on the geometry and grid 
resolution in the box. 



6.2 The influence of the lower impenetrable 
boundary 

The structure of the SAL in KC2 (or KC which is identi- 
cal) is at odds with the SAL constructed with the Canuto- 
Mazzitelli (1991) approach (see Fig 1. in Demarque, Guen- 
ther & Kim (1999)) or presented in the numerical simulation 
by Rosenthal et al. The most obvious difference is that KC2 
predicts a broader SAL with a maximum of about 0.4, while 
Canuto-Mazzitelli predicts a peak of about unity and Rosen- 
thal et al a peak of 0.6. The height of the peak is related 
to the convective efficiency with respect to radiative heat 
transport. 

Models KC2, A and B differ only in their domain 
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Table 1. List of Simulations 



Model 


Width (Mm) 


Depth (Mm) 


Depth (PSH) 


N x x N y x N z 


A* (km) 


A z (km) 


KC/KC2 


1.35 


0.9 


1.4 


60 3 


26 


17.5 


A 


1.35 


2.1 


7.2 


59 2 120 


26 


17.5 


B 


1.35 


2.8 


8.5 


59 2 170 


26 


17.5 


C 


2.7 


2.8 


8.5 


116 2 170 


26 


17.5 


D 


2.7 


2.8 


8.5 


58 2 170 


52 


17.5 


E 


5.4 


2.8 


8.5 


114 2 170 


52 


17.5 



depths. Simulation A was constructed by adding 3 extra 
pressure scale heights to the base of model KC2 (see Ta- 
ble 1). We first thermally relaxed the lower layer, because 
with only hydrostatic support (no turbulent pressure), the 
overlying layer collapsed, and the simulation crashed. Also, 
when joining the two layers, the turbulent viscosity was tem- 
porarily increased. This smoothed out the fluctuations that 
were produced by suddenly removing the lower boundary. 
Model B was made by adding a hydrostatic layer computed 
using the MLT to model A. This increased the depth from 
2.1 Mm to 2.8 Mm. As the turbulent pressure near the bot- 
tom of model A was small (about 1 % of the gas pressure) , 
in this case we did not need to relax the hydrostatic lower 
layer separately. 

6.2.1 Effect on kinetic energy and turbulent pressure 

To demonstrate the effect of the lower boundary on the tur- 
bulent flow, we computed the ratio of the horizontal kinetic 
energy to the gas pressure. This is equivalent to the hori- 
zontal turbulent Mach number squared, i.e. 

Xhoriz = -p{V x +V y )/P=-(v x +V y )/c a (22) 

Fig. |l4| shows Xhoriz for models KC2, A and B. In KC2 the 
magnitude of Xhoriz shoots up at the base. The net effect 
of the fast downflows striking the lower boundary was to 
speed up the overall flow. This increased the total turbulent 
kinetic energy x throughout the box (Fig. |l5| ). 

The upturn near the bottom of each layer, is clearly re- 
duced as the boundary is moved deeper, and is almost elimi- 
nated when the depth is 2.8 Mm. This problem is much less 
severe at the top boundary because the underlying region 
is sub-adiabatic (the radiative zone). The vertical velocity 
is already small when the flow hits the top boundary. The 
effect on Pturb of moving down the lower boundary is shown 
in Fig. [l(| In model KC2, P tur b drops sharply as the impen- 
etrable bottom is approached. While for A and B it has a 
much smoother decay. 

6.2.2 Effect on superadiabaticity 

The run of V - V ad for the MLT, and models KC2, and 
C was shown previously in Fig. ^ When the SAL is moved 
outwards by turbulent pressure, the convective efficiency is 
reduced and radiation is forced to carry more of the total 
flux. This should result in an increased height of the SAL 
compared to the MLT. However, in KC2 the convection is 
overly turbulent and is thus (incorrectly) more efficient than 
convection computed via the MLT. This results in a drop in 
the height of the SAL compared to the MLT. 



6.3 The influence of the upper impenetrable 
boundary 

As the top boundary is also impenetrable we need to ensure 
that its position does not effect the convection either. We 
need to prove that the decrease of vertical velocity between 
the convection surface (where the depth = Mm) and the 
top of the box, is not due to the top boundary. Rather it 
should be due to the stable layer at the top, i.e. convection- 
radiation losses. 

To address this issue we damped the horizontal veloc- 
ity at the top by replacing the stress free boundary with 
a no-slip top. The flow is then relaxed and statistics are 
gathered as usual. Fig. [j^ shows V — X7 a d,v" and v" for 
simulations with slip and no-slip top boundaries. The veloc- 
ities are non-dimensionalised by the local isothermal sound 
speed. This conveniently enables them to be on the same 
axis as V — V a d- Apart from the top 100 km, where the 
two w"'s diverge because of the different top boundary con- 
ditions, the rms vertical and horizontal velocities are nearly 
the same for the no-slip and stress free top boundaries. This 
implies the drop in vertical velocity near the top is primarily 
due to to convection to radiation losses, rather than the top 
boundary. Furthermore, the horizontal velocity in the upper 
atmosphere does not affect the SAL structure (i.e. V — V al j) 
much. 

6.4 The influence of the domain width 

Models B, C, D and E have widths of 1.35 Mm, 2.7 Mm, 2.7 
Mm and 5.4 Mm. To judge the effect of width on simula- 
tions without changing the grid spacing, we should compare 
B to C and D to E. All four models have similar Pturb and 
V — V a( j (Figs. ^| and ^) . Increasing the width seems to 
have a minor effect on either Pturb or V — V a d- In general it 
is essential to resolve turbulent motions inside the SAL re- 
gion. As this is only about 250 km thick, the grid spacing in 
the SAL need to be very small. To show that all our simula- 
tions have resolved the SAL, we have plotted the individual 
vertical grid points on the figure. 

The variation of x with domain width is more interest- 
ing. This is shown in Fig. When the width is increased 
from 1.35 Mm to 2.7 Mm, x increases (especially near the 
top). However, when the width is increased from 2.7 Mm 
to 5.4 Mm there is only a very small change in x- A box 
width of 2.7 Mm seems to be a sufficient to resolve x f° r 
solar granules. 

To provide further evidence that 2.7 Mm is a large 
enough box width, we computed v",v' y ' and for mod- 
els D and E. Fig. Ell shows all three velocity components for 
both D and E. Doubling the width has only a minor effect 



10 Robinson et al. 



on any particular velocity component. The small differences 
in the deeper part, are due to insufficient convergence. As 
the domain has periodic lateral boundaries, eventually the 
horizontal velocities should all be the same. 

Fig. shows contour plots of the instantaneous vertical 
velocity, for a horizontal cross-section, for models of different 
widths. The uppermost frame shows that not even a single 
granule can fit in the box when the width is 1.35 Mm, while 
about 2 fit easily in to the box in the second frame. In the 
final frame about 9 or 10 larger cells fit in to the box. Due to 
computational restrictions the last figure is computed on a 
coarser mesh than the previous three frames so the granules 
are not clearly depicted. 
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6.5 The influence of numerical resolution 

While the main point of these tests is to show that the walls 
of the computational box can have a significant effect on 
granular convection, we can also partially address the effect 
of changing the horizontal grid spacing. Models C and D 
produced very similar P tur b and V — V a d, while \ differed 
only slightly. As \ depends on the horizontal as well as the 
vertical velocity, this probably reflects the sensitivity of the 
horizontal component of kinetic energy to horizontal grid 
resolution. This is particularly noticeable near the top where 
the flow is mostly horizontal. 
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Figure 1. The rms turbulent velocities in the horizontal and ver- Figure 4. Supcradiabaticity versus fractional radius. The crosses 

tical directions versus depth. The closeness of the two horizontal are from the ID stellar model (MLT), the dashes are for model 

velocities confirms that the simulation is close to statistical con- KC2 and the triple dot-dash is for model C (see appendix for dc- 

vergence. tails). In both KC2 and C the original (MLT) convective bound- 

ary is moved out by turbulent pressure. 
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Figure 2. C[v' x v' y ] for three integration times. As C[v' x v' y ] con- 
verges more slowly than most other quantities, the degree of sta- 
tistical convergence can be estimated by how close C[v' x v' y ] is to 
zero. Convergence is slowest near the bottom of the box. 
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Figure 3. The mean and rms vertical velocity in model C. 
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Figure 5. The difference between the observed and the computed 
p-modc frequencies for a Standard Solar model, SSM (see text). 
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Figure 6. The difference between the observed and the computed 
p-modc frequencies for a solar model with turbulence included 
from model KC2 (see appendix). 
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Figure 7. The difference between the observed and the computed 
p-modc frequencies with turbulence included from model C (see 
appendix) . 




Figure 8. Contours of vertical velocity at one instant in time. 
The darker regions are downfiows, and the lighter regions are up- 
flows. From top to bottom, the frames are at depths of -0.14 Mm, 
0.03 Mm, 0.2 Mm and 1.0 Mm. Positive depths are measured 
inwards from the solar surface. The width of each frame is 3.75 
Mm. The small parallel vertical lines in the first frame indicate 
significant grid oscillations. By slightly increasing the SGS viscos- 
ity at the top, the oscillations are damped out before starting the 
actual statistical computations. The contours themselves suggest 
that the granules remain correlated throughout the SAL (which 
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Figure 9. The half width of the 2-point vertical correlation Figure 10. 2-point vertical velocity correlation at 5 different 

length of vertical velocity at different depths. The right side of depths. Unlike the case of deep convection, the velocity corrc- 

thc plateau is 3 PSH above the bottom of the box. Between the lation does not scale with pressure scale height 
surface and a depth of 1 Mm, the half-width is nearly constant. 
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Figure 11. 2-point entropy correlation at the same levels as in 
the previous figure. Close to the surface, ascending fluid parcels 
loose their entropy identity much sooner than descending parcels. 




Figure 12. Correlation coefficient between entropy and vertical 
velocity for upflows and downflows. Also plotted is V — V a( j for 
comparison. 
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Figure 13. Hydrostatic structure of two ID stellar models (see 
text) for KC (dotted) and KC2 (solid). The mean thermal struc- 
ture from the hydrodynamic simulations is also shown for KC (di- 
amonds) and KC2 (crosses). Despite different initial stellar mod- 
els, the turbulence causes both models to eventually converge to 
the same equilibrium state. 



Figure 16. Ratio of turbulent pressure to gas pressure, in this 
case denoted by Pturbi f° r simulations of different depths. 
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Figure 14. Non-dimensional horizontal turbulent kinetic energy 
per unit mass, for simulations of different depths. Note the speed- 
ing up of the flow at the lower stress free boundary in KC2. 
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Figure 17. Model D with a slip (stress free) and a no-slip top 
boundary. The figure includes horizontal and vertical turbulent 
velocities and superadiabaticity. The velocities have been scaled 
by the local sound speed. 




Figure 15. Non-dimensional turbulent kinetic energy per unit 
mass, x f° r simulations of different depths. 
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Figure 18. Turbulent pressure as a function of domain width, 
fturb is the ratio of turbulent pressure to gas pressure. The plot 
for D and E has been shifted by 0.5 Mm for clarity. Without the 
shift the peaks would coincide. 



Figure 21. The three components of turbulent velocity (scaled 
by the isothermal sound speed), for simulations D and E. As the 
change in width has little effect on any component of the velocity, 
2.7 Mm seems to be a sufficient box width. 
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Figure 19. Superadiabaticity as a function of domain width for 
models B, C, D and E. The vetical grid points are individually 
marked to show that the SAL is well resolved. The plots have 
been spaced by 0.5 for clarity. Without the horizontal spacing all 
the plots would have zero superadiabaticity at a depth of Mm. 
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Figure 20. Non-dimensional turbulent kinetic energy for differ- 
ent domain widths. The plot for D and E has been shifted by 0.5 
Mm for clarity. 
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Figure 22. Snapshot of vertical velocity contours for different 
domain widths. From top to bottom the domain widths are 1.35, 
2.7, 3.75 and 5.4 Mm, respectively. Due to computational restric- 
tions the last frame has half the resolution than the other 3. 



